####### FIGURE A1 ##########

reds.loc = read.csv("data/reds_locations_geocoded.csv", stringsAsFactors = F) %>%
  mutate_at(vars(contains("_code")), .funs = as.numeric)
dn.loc = read.csv("data/dunning_nilekani_locations_geocoded.csv", stringsAsFactors = F)

# Get India map
india.st = read_sf("data/india_states/STATE2011.shp")

map.data = reds %>%
  group_by(stateid, districtid, blockid, villageid) %>%
  summarise_at(vars(vwreserv_ever, vstreserv_ever, vscreserv_ever),
               .funs = max, na.rm = T) %>%
  left_join(reds.loc %>%
              dplyr::select(state_code, dist_code, block_code,village_code, lat, lon),
            by = c("stateid" = "state_code", "districtid" = "dist_code",
                   "blockid" = "block_code", "villageid" = "village_code")) %>%
  rename(reserved_fem = vwreserv_ever, reserved_sc = vscreserv_ever,
         reserved_st = vstreserv_ever) %>%
  mutate(dataset = "All India") %>%
  bind_rows(dn.dta.vill %>%
              dplyr::select(District_code, Block_code,  Panchayat_code, reserved_fem,
                            reserved_sc, reserved_st) %>%
              left_join(dn.loc %>% dplyr::select(District_code, Block_code, Panchayat_code, lat, lon),
                        by = c("District_code" = "District_code", "Block_code" = "Block_code",
                               "Panchayat_code" = "Panchayat_code")) %>%
              rename(districtid = District_code, blockid = Block_code,
                     villageid = Panchayat_code) %>%
              mutate(dataset = "BR-RJ-KA"))

# Create variable for states/UTs that are missing from the data
nodata = c("Andaman & Nicobar Islands", "Arunachal Pradesh", "Assam",
           "Chandigarh", "Dadra & Nagar Haveli", "Daman & Diu",
           "Delhi", "Goa", "Jammu & Kashmir", "Lakshadweep", "Manipur",
           "Meghalaya", "Mizoram", "Nagaland", "Puducherry",
           "Sikkim", "Tripura", "Uttarakhand")
india.st = india.st %>% mutate(has_data = ifelse(NAME %in% nodata, 0, 1))


india = st_as_sf(india.st)

# Female reservations
pl1 = ggplot()+
  geom_sf(data = india, aes(fill = as.factor(has_data)), color = "black", lwd = 0.2, alpha = 0.5)+
  geom_point(data = map.data %>% filter(!(is.na(reserved_fem)|is.infinite(reserved_fem))&dataset=="All India"),
             aes(x = lon, y = lat, shape = as.factor(reserved_fem), color = as.factor(reserved_fem)),
             size = 2)+
  scale_shape_manual(values = c(17, 16), name = "Women's Quota",
                     labels = c("No", "Yes"))+
  scale_color_manual(values = c("turquoise", "firebrick4"), name = "Women's Quota",
                     labels = c("No", "Yes"))+
  scale_fill_manual(values = c("gray32", "white"), name = "State With REDS Villages",
                    labels = c("No", "Yes"), guide = "none")+
  theme_pubr()+
  theme(axis.line=element_blank(),axis.text.x=element_blank(),
        axis.text.y=element_blank(),axis.ticks=element_blank(),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),legend.position="top",
        panel.background=element_blank(),panel.border=element_blank(),panel.grid.major=element_blank(),
        panel.grid.minor=element_blank(),plot.background=element_blank(), 
        legend.title = element_text(size = 20), 
        legend.text = element_text(size = 20))

pl1  
ggsave(filename = paste0(fig.out, "figure_A1a.png"))

# SC reservations
pl2 = ggplot(data = india)+
  geom_sf(aes(fill = as.factor(has_data)), color = "black", lwd = 0.2, alpha = 0.5)+
  geom_point(data = map.data %>% filter(!(is.na(reserved_sc)|is.infinite(reserved_sc))&dataset=="All India"),
             aes(x = lon, y = lat, shape = as.factor(reserved_sc), color = as.factor(reserved_sc)),
             size = 2)+
  scale_shape_manual(values = c(17, 16), name = "SC Quota",
                     labels = c("No", "Yes"))+
  scale_color_manual(values = c("turquoise", "firebrick4"), name = "SC Quota",
                     labels = c("No", "Yes"))+
  scale_fill_manual(values = c("gray32", "white"), name = "State With REDS Villages",
                    labels = c("No", "Yes"), guide = "none")+
  theme_pubr()+
  theme(axis.line=element_blank(),axis.text.x=element_blank(),
        axis.text.y=element_blank(),axis.ticks=element_blank(),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),legend.position="top",
        panel.background=element_blank(),panel.border=element_blank(),panel.grid.major=element_blank(),
        panel.grid.minor=element_blank(),plot.background=element_blank(),
        legend.title = element_text(size = 20), 
        legend.text = element_text(size = 20))
pl2
ggsave(filename = paste0(fig.out, "figure_A1c.png"))

# ST reservations
pl3 = ggplot(data = india)+
  geom_sf(aes(fill = as.factor(has_data)), color = "black", lwd = 0.2, alpha = 0.5)+
  geom_point(data = map.data %>% filter(!(is.na(reserved_st)|is.infinite(reserved_st))&dataset=="All India"),
             aes(x = lon, y = lat, shape = as.factor(reserved_st),
                 color = as.factor(reserved_st)),
             size = 2)+
  scale_shape_manual(values = c(17, 16), name = "ST Quota",
                     labels = c("No", "Yes"))+
  scale_color_manual(values = c("turquoise", "firebrick4"), name = "ST Quota",
                     labels = c("No", "Yes"))+
  scale_fill_manual(values = c("gray32", "white"), name = "State With REDS Villages",
                    labels = c("No", "Yes"), guide = "none")+
  theme_pubr()+
  theme(axis.line=element_blank(),axis.text.x=element_blank(),
        axis.text.y=element_blank(),axis.ticks=element_blank(),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),legend.position="top",
        panel.background=element_blank(),panel.border=element_blank(),panel.grid.major=element_blank(),
        panel.grid.minor=element_blank(),plot.background=element_blank(),
        legend.title = element_text(size = 20), 
        legend.text = element_text(size = 20))
pl3
ggsave(filename = paste0(fig.out, "figure_A1b.png"))



############ FIGURE A2 ###############
reds.vill = reds %>%
  group_by(stateid, districtid, tehsilid,
           blockid, villageid, period_new) %>%
  summarise_at(vars(vwreserv, vscreserv, vstreserv, vscstreserv,
                    veducilo_5670, cult_prop_99, distr_HQ,
                    prop_SC91, prop_ST91, vpop_91, vpest_99,
                    vflood_99, vdrought_99), .funs = mean, na.rm = T) %>%
  mutate(vpop_91_10K = ifelse(is.na(vpop_91), NA, vpop_91/10000),
         distr_HQ = ifelse(is.na(distr_HQ), NA, distr_HQ/100),
         vscstreserv = ifelse(vscreserv==1 | vstreserv == 1, 1, 0))

reds.vill = reds.vill %>%
  group_by(stateid, districtid, tehsilid, blockid, villageid) %>%
  mutate(no_res = ifelse(vwreserv==0 & vstreserv==0, 1, 0),
         vwreserv_only = ifelse(vwreserv==1 & vstreserv==0, 1, 0),
         vstreserv_only = ifelse(vwreserv==0 & vstreserv==1, 1, 0),
         vwreserv_andst = ifelse(vwreserv==1&vstreserv==1, 1, 0),
         no_res_sum = sum(no_res, na.rm = T),
         vwreserv_only_sum = sum(vwreserv_only, na.rm = T),
         vstreserv_only_sum = sum(vstreserv_only, na.rm = T),
         vwreserv_andst_sum = sum(vwreserv_andst, na.rm = T))

treatment = reds.vill %>%
  filter(!(is.na(vstreserv)|is.na(vwreserv))) %>%
  ungroup() %>%
  mutate(reservation = case_when(
    vstreserv==0 & vwreserv==0 ~ 1,
    vstreserv==0 & vwreserv==1 ~ 2,
    vstreserv==1 & vwreserv==0 ~ 3,
    vstreserv==1 & vwreserv==1 ~ 4
  ),
  reservation_fac = factor(reservation, labels = c("No quota",
                                                   "Women's quota",
                                                   "ST quota",
                                                   "Women's & ST quota")),
  prop_ST91_bin = case_when(
    prop_ST91 < 0.1 ~ 0.0,
    prop_ST91 >=0.1 & prop_ST91 < 0.2 ~ 0.1,
    prop_ST91 >=0.2 & prop_ST91 < 0.3 ~ 0.2,
    prop_ST91 >=0.3 & prop_ST91 < 0.4 ~ 0.3,
    prop_ST91 >=0.4 & prop_ST91 < 0.5 ~ 0.4,
    prop_ST91 >=0.5 & prop_ST91 < 0.6 ~ 0.5,
    prop_ST91 >=0.6 & prop_ST91 < 0.7 ~ 0.6,
    prop_ST91 >=0.7 & prop_ST91 < 0.8 ~ 0.7
  )) %>%
  group_by(prop_ST91_bin, reservation_fac) %>%
  add_count(name = "count") %>%
  ungroup %>%
  group_by(prop_ST91_bin) %>%
  add_count(name = "reservation_count") %>%
  mutate(proportion = count/reservation_count)%>%
  group_by(reservation_fac, prop_ST91_bin) %>%
  summarise_at(vars(proportion), .funs = mean) %>%
  mutate(label = round(proportion, 2) %>% as.character()) %>%
  group_by(prop_ST91_bin) %>%
  arrange(prop_ST91_bin, desc(reservation_fac)) %>%
  mutate(cumsum = cumsum(proportion)) %>%
  mutate(lag_cumsum = dplyr::lag(cumsum, 1)) %>%
  mutate(lag_cumsum = ifelse(is.na(lag_cumsum), 0, lag_cumsum)) %>%
  mutate(pos = (cumsum+lag_cumsum)/2) %>%
  mutate(label = ifelse(proportion < 0.05, "", label))
treat = ggplot(data = treatment, aes(x = as.factor(prop_ST91_bin), y = proportion, label = label))+
  geom_bar(aes(fill = reservation_fac), stat = "identity", position = "fill")+
  geom_text(aes(y = pos, color = ifelse(reservation_fac %in% c("Not reserved", "Women's reservation"), "black", "white")), size = 4)+
  scale_fill_manual(values = c("gray92", "gray75", "gray53", "firebrick4"),
                    name = "Quota type")+
  scale_y_continuous(labels=scales::percent, limits = c(0, 1))+
  scale_x_discrete(labels = c("0-9%", "10-19%", "20-29%", "30-39%", "40-49%",
                              "50-59%", "60-69%"))+
  scale_color_manual(values = c("black", "white"), guide = "none")+
  labs(x = "Proportion of Scheduled Tribes (1991)",
       y = "Proportion of villages")+
  theme_pubr()+
  theme(axis.text = element_text(face = "bold", size = 16),
        axis.title = element_text(face = "bold", size = 17),
        legend.title = element_text(face = "bold", size = 16),
        legend.text = element_text(face = "bold", size = 15))
treat
ggsave(treat, filename = paste0(fig.out, "figure_A2.png"))




###### FIGURE A4: BALANCE TABLES #######
labels = c( "% Cultivators", "% Female Workers",
            "% Female",
            "% Without education", "% Literate", "% Marginal workers",
            "% Scheduled Caste", "% Scheduled Tribe",
            "% Working population")


###### SC/ST VS. NONE
dn.dta$reserved_femscst = ifelse(dn.dta$reserved_fem==1 & (dn.dta$reserved_sc==1|dn.dta$reserved_st==1), 1, 0)
dn.dta$reserved_scst = ifelse((dn.dta$reserved_sc==1|dn.dta$reserved_st==1), 1, 0)
dn.dta.vill$reserved_femscst = ifelse(dn.dta.vill$reserved_fem==1 & (dn.dta.vill$reserved_sc==1|dn.dta.vill$reserved_st==1), 1, 0)
dn.dta.vill$reserved_scst = ifelse((dn.dta.vill$reserved_sc==1|dn.dta.vill$reserved_st==1), 1, 0)

bal.scst = bal_tab(data = dn.dta.vill,
                   group.var = reserved_scst,
                   balancevar = pc_sc, pc_st, pc_fem,
                   pc_lit, pc_ill, pc_work, pc_f_work, pc_margwork, pc_cult)

bal.fem = bal_tab(data = dn.dta.vill,
                  group.var = reserved_fem,
                  balancevar = pc_sc, pc_st, pc_fem,
                  pc_lit, pc_ill, pc_work, pc_f_work, pc_margwork, pc_cult)


bal.dta = bind_rows(bal.scst[[1]] %>% mutate(sample = "SC/ST vs. No Quota"),
                    bal.fem[[1]] %>% mutate(sample = "Women's vs. No Quota"))

fig = bal_plot(data = bal.dta, labs = labels)+
  theme_pubr()+
  theme(axis.title = element_text(face = "bold", size = 20),
        axis.text = element_text(face = "bold", size = 18),
        legend.title = element_text(face = "bold", size = 18),
        strip.text = element_text(face = "bold", size = 18),
        legend.text = element_text(size = 18))
fig
ggsave(fig, filename = paste0(fig.out, "figure_A4.png"), 
       width = 12.8, height = 7.5)


lm(reserved_fem ~ pc_sc+pc_st+pc_fem+pc_lit+pc_ill+pc_work+pc_f_work+pc_margwork+pc_cult,
   data = dn.dta.vill) %>% summary()
# F-statistic: 0.7159; p-value: 0.6774

lm(reserved_scst ~ pc_sc+pc_st+pc_fem+pc_lit+pc_ill+pc_work+pc_f_work+pc_margwork+pc_cult,
   data = dn.dta.vill) %>% summary()
# F-statistic: 0.6551; p-value: 0.7493



####### FIGURE A.3: REDS MATCHED - ST RESERVATIONS #######
reds.st.matched.vill = reds.st.matched %>%
  group_by(stateid, districtid, tehsilid,
           blockid, villageid, period_new) %>%
  summarise_at(vars(vwreserv, vscreserv, vstreserv,
                    veducilo_5670, cult_prop_99, distr_HQ,
                    prop_SC91, prop_ST91, vpop_91, vpest_99,
                    vflood_99, vdrought_99), .funs = mean, na.rm = T) %>%
  mutate(vpop_91_10K = ifelse(is.na(vpop_91), NA, vpop_91/10000),
         distr_HQ = ifelse(is.na(distr_HQ), NA, distr_HQ/100),
         vscstreserv = ifelse(vscreserv==1 | vstreserv == 1, 1, 0))

reds.scst.matched.vill = reds.scst.matched %>%
  group_by(stateid, districtid, tehsilid,
           blockid, villageid, period_new) %>%
  summarise_at(vars(vwreserv, vscreserv, vstreserv,
                    veducilo_5670, cult_prop_99, distr_HQ,
                    prop_SC91, prop_ST91, vpop_91, vpest_99,
                    vflood_99, vdrought_99), .funs = mean, na.rm = T) %>%
  mutate(vpop_91_10K = ifelse(is.na(vpop_91), NA, vpop_91/10000),
         distr_HQ = ifelse(is.na(distr_HQ), NA, distr_HQ/100),
         vscstreserv = ifelse(vscreserv==1 | vstreserv == 1, 1, 0))

labels = c("% Cultivated area in 1999", "Distance to district HQ (100Km)",
           "% SC (1991)", "% ST (1991)", "Any drought in 1999", "% No or incomplete education",
           "Any flood in 1999", "Any pest in 1999")



###### SC/ST & NONE
bal.scst = bal_tab(data = reds.vill %>% ungroup(),
                   group.var = vscstreserv,
                   balancevar = veducilo_5670, cult_prop_99, distr_HQ,
                   prop_SC91, prop_ST91, vpest_99,
                   vflood_99, vdrought_99)

bal.scst.scstmatched = bal_tab(data = reds.scst.matched.vill %>% ungroup(),
                               group.var = vscstreserv,
                               balancevar = veducilo_5670, cult_prop_99, distr_HQ,
                               prop_SC91, prop_ST91, vpest_99,
                               vflood_99, vdrought_99)

bal.scst.stmatched = bal_tab(data = reds.st.matched.vill %>% ungroup(),
                             group.var = vscstreserv,
                             balancevar = veducilo_5670, cult_prop_99, distr_HQ,
                             prop_SC91, prop_ST91, vpest_99,
                             vflood_99, vdrought_99)



bal.dta1 = bind_rows(
  bal.scst[[1]] %>% mutate(sample = "REDS"),
  bal.scst.stmatched[[1]] %>% mutate(sample = "REDS (Matched: ST)"),
  bal.scst.scstmatched[[1]] %>% mutate(sample = "REDS (Matched: SC/ST)")
)

fig = bal_plot(data = bal.dta1, labs = labels)+
  theme_pubr()+
  theme(axis.title = element_text(face = "bold", size = 20),
        axis.text = element_text(face = "bold", size = 18),
        legend.title = element_text(face = "bold", size = 18),
        strip.text = element_text(face = "bold", size = 16),
        legend.text = element_text(size = 18))
fig
ggsave(fig, filename = paste0(fig.out, "figue_A3a.png"), width = 17.7,
       height = 7.5)





###### FEM VS. NONE
bal.fem = bal_tab(data = reds.vill %>% ungroup(),
                  group.var = vwreserv,
                  balancevar = veducilo_5670, cult_prop_99, distr_HQ,
                  prop_SC91, prop_ST91, vpest_99,
                  vflood_99, vdrought_99)

bal.fem.scmatched = bal_tab(data = reds.scst.matched.vill %>% ungroup(),
                            group.var = vwreserv,
                            balancevar = veducilo_5670, cult_prop_99, distr_HQ,
                            prop_SC91, prop_ST91, vpest_99,
                            vflood_99, vdrought_99)

bal.fem.stmatched = bal_tab(data = reds.st.matched.vill %>% ungroup(),
                            group.var = vwreserv,
                            balancevar = veducilo_5670, cult_prop_99, distr_HQ,
                            prop_SC91, prop_ST91, vpest_99,
                            vflood_99, vdrought_99)



bal.dta1 = bind_rows(
  bal.fem[[1]] %>% mutate(sample = "REDS"),
  bal.fem.scmatched[[1]] %>% mutate(sample = "REDS (Matched: SC/ST)"),
  bal.fem.stmatched[[1]] %>% mutate(sample = "REDS (Matched: ST)")
)
fig = bal_plot(data = bal.dta1, labs = labels)+
  theme_pubr()+
  theme(axis.title = element_text(face = "bold", size = 20),
        axis.text = element_text(face = "bold", size = 18),
        legend.title = element_text(face = "bold", size = 18),
        strip.text = element_text(face = "bold", size = 16),
        legend.text = element_text(size = 18))
fig


ggsave(fig, filename = paste0(fig.out, "figure_A3b.png"), width = 17.7,
       height = 7.5)






########### TABLE B.3 ################


#### DUNNING AND NILEKANI
varnames = c("marry_other_caste_bin", "marry_caste_kar", "marry_other_caste_kar",
             "marry_caste_br", "marry_other_caste_br", "age", "vote_panch",
             "reserved_fem",  "reserved_st", "reserved_scst", "reserved_femst",
             "reserved_femscst", "prop_st")
st_fem_sum = summary.fun(data = dn.dta %>% filter(fem==1 & st==1), varnames = varnames)
all_fem_sum = summary.fun(data = dn.dta %>% filter(fem==1), varnames = varnames)
entire_sum = summary.fun(data = dn.dta, varnames = varnames)


labels = c("Marry Other Caste (0/1)",
           "Marry Other Caste - KN (1-5)", "Marry Other Caste - KN (0/1)",
           "Marry Other Caste - Raj-Bih (1-7)", "Marry Other Caste - Raj-Bih (0/1)",
           "Age",
           "Voted in last panchayat election (0/1)", "Women's Reservations (0/1)",
           "ST Reservations (0/1)", "ST/SC Reservations (0/1)", "Women's and ST Reservations (0/1)",
           "Women's and ST/SC Reservations (0/1)", 
           "% ST in village")
descr = entire_sum %>% dplyr::select(mean, sd) %>%
  rename(mean_ent = mean, sd_ent = sd) %>%
  bind_cols(st_fem_sum %>% dplyr::select(mean, sd) %>%
              rename(mean_stf = mean, sd_stf = sd)) %>%
  bind_cols(all_fem_sum %>% dplyr::select(mean, sd) %>%
              rename(mean_allf = mean, sd_allf = sd)) %>%
  mutate_all(.funs = round, digits = 3) %>%
  mutate(var = labels)  %>%
  dplyr::select(var, contains("ent"), contains("stf"), contains("allf"))


print(xtable(descr, digits = 2), only.contents = T, include.colnames = F,
      booktabs = F, include.rownames = F, hline.after = NULL,
      file = paste0(tab.out, "tableB3_panelB.tex"))


### REDS
reds = reds %>%
  mutate(caste_interact_not_easy = ifelse(caste_interact_sc>=0.2 & caste_interact_sc<0.8, 1, 0)) %>%
  mutate(caste_interact_not_easy = ifelse(is.na(caste_interact_sc), 1, caste_interact_not_easy)) %>%
  mutate(caste_interact_easy_2 = ifelse(caste_interact_sc==0.8, 1, 0)) %>%
  mutate(caste_interact_easy_2 = ifelse(is.na(caste_interact_easy_2), 0, caste_interact_easy_2)) %>%
  mutate(caste_interact_very_easy = ifelse(caste_interact_sc==1, 1, 0)) %>%
  mutate(caste_interact_very_easy = ifelse(is.na(caste_interact_very_easy), 0, caste_interact_very_easy)) %>%
  mutate(vwreservXvstreserv = ifelse(vwreserv==1 | vstreserv==1, 1, 0)) %>%
  mutate(vwreservXvscstreserv = ifelse(vwreserv==1 | vscstreserv==1, 1, 0))
varnames = c("caste_interact_bi_comp", "caste_interact_bi", "caste_interact_not_easy",
              "caste_interact_easy_2",
             "caste_interact_very_easy", "top20land",
             "birthyr", "fem", "Caste_ST", "Caste_SC",
             "prop_STnow", "prop_SCnow",
             "vwreserv", "vstreserv", "vscstreserv",
             "vwreservXvstreserv", "vwreservXvscstreserv",
             "caste_based_inter_family_disputes",
             "prog_indben", "prop_ST91", "prop_SC91",
             "SC_party_pc_vids")
reds_st_fem_sum = summary.fun(data = reds %>% filter(fem==1 & Caste_ST==1), varnames = varnames)
reds_all_fem_sum = summary.fun(data = reds %>% filter(fem==1), varnames = varnames)
reds_entire_sum = summary.fun(data = reds, varnames = varnames)

labels = c("Easy to interact with other castes (0/1)",
           "Easy to interact (0/1) - missing dropped",
           "Interacting with other castes: Not easy (0/1)",
           "Interacting with other castes: Easy (0/1)",
           "Interacting with other castes: Very easy (0/1)",
           "Top 20% landholder (0/1)", "Birth year", "Female (0/1)",
           "Caste: ST", "Caste: SC",
           "% ST in village today", "% SC in village today",
           "Women's Reservations (0/1)",
           "ST Reservations (0/1)", "SC/ST Reservations (0/1)", "Women's and ST Reservations (0/1)",
           "Women's and SC/ST Reservations (0/1)",
           "Inter-family Caste-Based Conflict",
           "Receive benefits (0/1)",
           "% ST in subdistrict (1991)",
           "% SC in subdistrict (1991)", 
           "ST Party Vote Share in Vidhan Sabha")

descr = reds_entire_sum %>% dplyr::select(mean, sd) %>%
  rename(mean_ent = mean, sd_ent = sd) %>%
  bind_cols(reds_st_fem_sum %>% dplyr::select(mean, sd) %>%
              rename(mean_stf = mean, sd_stf = sd)) %>%
  bind_cols(reds_all_fem_sum %>% dplyr::select(mean, sd) %>%
              rename(mean_allf = mean, sd_allf = sd)) %>%
  mutate_all(.funs = round, digits = 3) %>%
  mutate(var = labels) %>%
  dplyr::select(var, contains("ent"), contains("stf"), contains("allf"))

print(xtable(descr, digits = 2), only.contents = T, include.colnames = F,
      booktabs = F, include.rownames = F, hline.after = NULL,
      file = paste0(tab.out, "TableB3_panelA.tex"))



######### FIGURE D6 ########

reds.pradhan = read_dta("data/vdeck36.dta") %>%
  mutate(districtid = districtid %>% destring, 
         tehsilid = tehsilid %>% destring,
         blockid = blockid %>% destring, 
         villageid = villageid %>% destring, 
         period_new = 4-q1) %>%
  group_by(districtid, tehsilid, blockid, villageid, period_new) %>%
  mutate(candidate_rank = order(desc(q12))) %>%
  filter(candidate_rank==1) %>%
  left_join(reds %>%
              group_by(districtid, tehsilid, blockid, villageid, period_new) %>%
              summarise_at(vars(vwreserv, vstreserv, vscreserv, vobcreserv), .funs = mean, na.rm = T) %>%
              mutate(reservation_type = case_when(
                vstreserv==0 & vwreserv==0  ~ "No reservation",
                vstreserv==0 & vwreserv==1 ~ "Women only",
                vstreserv==1 & vwreserv==0 ~ "ST only",
                vstreserv==1 & vwreserv==1 ~ "Women & ST"))) %>%
  mutate(woman = ifelse(q5==2, 1, 0), 
         man = ifelse(q5==1, 1, 0), 
         ST = ifelse(q8==2, 1, 0), 
         OC = ifelse(q8==4, 1, 0), 
         SC = ifelse(q8==1, 1, 0), 
         OBC = ifelse(q8==3, 1, 0), 
         ST_SC = ifelse(ST==1|SC==1, 1, 0)) %>%
  group_by(reservation_type) %>%
  summarise_at(vars(woman, man, ST, OC, SC, OBC, ST_SC), .funs = mean, na.rm = T)

reds.pradhan.sc = read_dta("data/vdeck36.dta") %>%
  mutate(districtid = districtid %>% destring, 
         tehsilid = tehsilid %>% destring,
         blockid = blockid %>% destring, 
         villageid = villageid %>% destring, 
         period_new = 4-q1) %>%
  group_by(districtid, tehsilid, blockid, villageid, period_new) %>%
  mutate(candidate_rank = order(desc(q12))) %>%
  filter(candidate_rank==1) %>%
  left_join(reds %>%
              group_by(districtid, tehsilid, blockid, villageid, period_new) %>%
              summarise_at(vars(vwreserv, vstreserv, vscreserv, vobcreserv), .funs = mean, na.rm = T) %>%
              mutate(reservation_type = case_when(
                vscreserv==0 & vwreserv==0  ~ "No reservation",
                vscreserv==0 & vwreserv==1 ~ "Women only",
                vscreserv==1 & vwreserv==0 ~ "SC only",
                vscreserv==1 & vwreserv==1 ~ "Women & SC",
                TRUE ~ "Other"))) %>%
  mutate(woman = ifelse(q5==2, 1, 0), 
         man = ifelse(q5==1, 1, 0), 
         ST = ifelse(q8==2, 1, 0), 
         OC = ifelse(q8==4, 1, 0), 
         SC = ifelse(q8==1, 1, 0), 
         OBC = ifelse(q8==3, 1, 0)) %>%
  group_by(reservation_type) %>%
  summarise_at(vars(woman, man, ST, OC, SC, OBC), .funs = mean, na.rm = T)


reds.pradhan.scst = read_dta("data/vdeck36.dta") %>%
  mutate(districtid = districtid %>% destring, 
         tehsilid = tehsilid %>% destring,
         blockid = blockid %>% destring, 
         villageid = villageid %>% destring, 
         period_new = 4-q1) %>%
  group_by(districtid, tehsilid, blockid, villageid, period_new) %>%
  mutate(candidate_rank = order(desc(q12))) %>%
  filter(candidate_rank==1) %>%
  left_join(reds %>%
              group_by(districtid, tehsilid, blockid, villageid, period_new) %>%
              summarise_at(vars(vwreserv, vscstreserv), .funs = mean, na.rm = T) %>%
              mutate(reservation_type = case_when(
                vscstreserv==0 & vwreserv==0  ~ "No reservation",
                vscstreserv==0 & vwreserv==1 ~ "Women only",
                vscstreserv==1 & vwreserv==0 ~ "SC/ST only",
                vscstreserv==1 & vwreserv==1 ~ "Women & SC/ST"))) %>%
  mutate(woman = ifelse(q5==2, 1, 0), 
         man = ifelse(q5==1, 1, 0), 
         ST_SC = ifelse(q8==2|q8==1, 1, 0), 
         OC = ifelse(q8==4, 1, 0), 
         OBC = ifelse(q8==3, 1, 0)) %>%
  group_by(reservation_type) %>%
  summarise_at(vars(woman, man, ST_SC, OC, OBC), .funs = mean, na.rm = T)


# gender
gender.pl = reds.pradhan %>%
  reshape2::melt(id.vars = "reservation_type") %>%
  filter(variable=="man"|variable=="woman") %>%
  mutate(label = round(value, 2) %>% as.character()) %>%
  group_by(reservation_type) %>%
  arrange(reservation_type, desc(variable)) %>%
  mutate(cumsum = cumsum(value)) %>%
  mutate(lag_cumsum = dplyr::lag(cumsum, 1)) %>%
  mutate(lag_cumsum = ifelse(is.na(lag_cumsum), 0, lag_cumsum)) %>%
  mutate(pos = (cumsum+lag_cumsum)/2) %>%
  mutate(label = ifelse(value < 0.05, "", label)) %>%
  ggplot(aes(x = reservation_type, y = value, label = label), group = variable)+
  geom_bar(aes(fill = as.factor(variable)),
           stat = "identity")+
  geom_text(aes(y = pos, color = ifelse(variable %in% c("man"), "black", "white")), size = 4)+
  scale_fill_manual(values = c("firebrick4", "gray"), name = "Gender")+
  scale_y_continuous(labels = scales::percent_format(accuracy = 1))+
  scale_color_manual(values = c("black", "white"), guide = "none")+
  labs(x= "Quota Type", y = "Percent of Pradhans")+
  theme_pubr()+
  theme(axis.text = element_text(face = "bold", size = 16),
        axis.title = element_text(face = "bold", size = 17),
        legend.title = element_text(face = "bold", size = 16),
        legend.text = element_text(face = "bold", size = 15))
gender.pl
ggsave(gender.pl, filename = paste0(fig.out, "figure_D6a.png"))

gender.pl.scst = reds.pradhan.scst %>%
  reshape2::melt(id.vars = "reservation_type") %>%
  filter(variable=="man"|variable=="woman") %>%
  mutate(label = round(value, 2) %>% as.character()) %>%
  group_by(reservation_type) %>%
  arrange(reservation_type, desc(variable)) %>%
  mutate(cumsum = cumsum(value)) %>%
  mutate(lag_cumsum = dplyr::lag(cumsum, 1)) %>%
  mutate(lag_cumsum = ifelse(is.na(lag_cumsum), 0, lag_cumsum)) %>%
  mutate(pos = (cumsum+lag_cumsum)/2) %>%
  mutate(label = ifelse(value < 0.05, "", label)) %>%
  ggplot(aes(x = reservation_type, y = value, label = label), group = variable)+
  geom_bar(aes(fill = as.factor(variable)),
           stat = "identity")+
  geom_text(aes(y = pos, color = ifelse(variable %in% c("man"), "black", "white")), size = 4)+
  scale_fill_manual(values = c("firebrick4", "gray"), name = "Gender")+
  scale_y_continuous(labels = scales::percent_format(accuracy = 1))+
  scale_color_manual(values = c("black", "white"), guide = "none")+
  labs(x= "Quota Type", y = "Percent of Pradhans")+
  theme_pubr()+
  theme(axis.text = element_text(face = "bold", size = 16),
        axis.title = element_text(face = "bold", size = 17),
        legend.title = element_text(face = "bold", size = 16),
        legend.text = element_text(face = "bold", size = 15))
gender.pl.scst
ggsave(gender.pl.scst, filename = paste0(fig.out, "figure_D6c.png"))


# caste - ST compliance
caste.pl = reds.pradhan %>%
  reshape2::melt(id.vars = "reservation_type") %>%
  filter(variable!="ST_SC") %>%
  filter(variable!="man"&variable!="woman") %>%
  group_by(reservation_type) %>%
  arrange(reservation_type, desc(variable)) %>%
  mutate(label = round(value, 2) %>% as.character()) %>%
  mutate(cumsum = cumsum(value)) %>%
  mutate(lag_cumsum = dplyr::lag(cumsum, 1)) %>%
  mutate(lag_cumsum = ifelse(is.na(lag_cumsum), 0, lag_cumsum)) %>%
  mutate(pos = (cumsum+lag_cumsum)/2) %>%
  mutate(label = ifelse(value < 0.05, "", label)) %>%
  ggplot(aes(x = reservation_type, y = value, label = label), group = variable)+
  geom_bar(aes(fill = as.factor(variable)),
           stat = "identity")+
  geom_text(aes(y = pos, color = ifelse(variable %in% c("OBC", "SC"), "black", "white")), size = 4)+
  scale_fill_manual(values = c("firebrick4", "gray47", "gray60", "gray80"),
                    name = "Caste")+
  scale_y_continuous(labels = scales::percent_format(accuracy = 1))+
  scale_color_manual(values = c("black", "white"), guide = "none")+
  labs(x= "Quota Type", y = "Percent of Pradhans")+
  theme_pubr()+
  theme(axis.text = element_text(face = "bold", size = 16),
        axis.title = element_text(face = "bold", size = 17),
        legend.title = element_text(face = "bold", size = 16),
        legend.text = element_text(face = "bold", size = 15))
caste.pl
ggsave(caste.pl, filename = paste0(fig.out, "figure_D6b.png"))

# caste - SC compliance
caste.pl.scst = reds.pradhan.scst %>%
  reshape2::melt(id.vars = "reservation_type") %>%
  filter(variable!="man"&variable!="woman") %>%
  group_by(reservation_type) %>%
  arrange(reservation_type, desc(variable)) %>%
  mutate(label = round(value, 2) %>% as.character()) %>%
  mutate(cumsum = cumsum(value)) %>%
  mutate(lag_cumsum = dplyr::lag(cumsum, 1)) %>%
  mutate(lag_cumsum = ifelse(is.na(lag_cumsum), 0, lag_cumsum)) %>%
  mutate(pos = (cumsum+lag_cumsum)/2) %>%
  mutate(label = ifelse(value < 0.05, "", label)) %>%
  ggplot(aes(x = reservation_type, y = value, label = label), group = variable)+
  geom_bar(aes(fill = as.factor(variable)),
           stat = "identity")+
  geom_text(aes(y = pos, color = ifelse(variable %in% c("OBC", "SC"), "black", "white")), size = 4)+
  scale_fill_manual(values = c("firebrick4", "gray47", "gray60"),
                    name = "Caste")+
  scale_y_continuous(labels = scales::percent_format(accuracy = 1))+
  scale_color_manual(values = c("black", "white"), guide = "none")+
  labs(x= "Quota Type", y = "Percent of Pradhans")+
  theme_pubr()+
  theme(axis.text = element_text(face = "bold", size = 16),
        axis.title = element_text(face = "bold", size = 17),
        legend.title = element_text(face = "bold", size = 16),
        legend.text = element_text(face = "bold", size = 15))
caste.pl.scst
ggsave(caste.pl.scst, filename = paste0(fig.out, "figure_D6d.png"))








